%%%Main File for BCDR Simulations - gains from industrial policy and trade
%%%policy

%This file computes the gains from a given trade and industrial policy
%using the SOE assumption

%CES Version

%Uses alternative TE assumptions (Figure 4 in BCDR)

clear all;

%digits(64);

K=32; %number of sectors
N=62; %includes ROW
K_man=15;    %number of manufacturing sectors

%%%Load data
%%%Data is organized as follows: first column origin country codes, second column
%%%sector codes, third column sectoral value added, fourth column aggregate origin trade balance,
%%%fifth column origin-sector expenditurs, 6th column starts the bilateral
%%%trade flows

%Note: data is sorted by sector, then origin.

A=importdata('simulation_ge_data.csv');
B=importdata('alternative_te_parameters.csv'); %these are the estimated alternative trade and scale elasticities elasticities
C=importdata('rho.csv');    %These are the re-estimated upper tier EoS, using the corresponding alternative TE estimates

%Data vectors
X_ink = A(:,6:5+N); %Matrix of trade flows, (NxK)xN
D_i=A(1:N,4); %vector of trade balances, Nx1
V_ik=A(:,3); %vector of sectoral value added, (NxK)x1
expenditures_ik=A(:,5); %vector of expenditures, (NxK)x1
countryid=A(1:N,1); %vector of country ids, Nx1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Reshaping data vectors

%Reshape trade matrix to be NxNxK (3 dimensional array)
%rows are exporters, columns are importers, 3rd dimension is sectors
X=zeros(N,N,K);
for k=1:K
    X(:,:,k) = X_ink(1+(k-1)*N:k*N,:);
end
X_ink=X;

%Other data matrices
V_ik = reshape(V_ik, [N 1 K]);
V_i=sum(V_ik,3);
D_i =reshape(D_i, [N 1 1]);
expenditures_ik = reshape(expenditures_ik,[N,1,K]);
expenditures_nk=permute(expenditures_ik,[2 1 3]);

%Create trade shares, expenditure shares
expenditures_nnk=repmat(expenditures_nk,[N 1 1]);
lambda_ink = X_ink./expenditures_nnk; %trade shares

expenditures_n=sum(expenditures_nk,3);
beta_nk=expenditures_nk./repmat(expenditures_n,[1 1 K]);

%Create sectoral domestic trade shares
dom_lambda_ik=zeros(N,K);
for i=1:N
    for k=1:k
        dom_lambda_ik(i,k) = lambda_ink(i,i,k);
    end
end
        
%Create total sectoral exports
exports_ik=zeros(N,K);
for i=1:N
    for k=1:k
        exports_ik(i,k) = V_ik(i,k)- X_ink(i,i,k);
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Gains from industrial policy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Het TE, Het SE (service 0 SE)
counter=0;
for j=1:11
    %Parameter vectors
    theta_nm=4.32; %non-manufacturing sector TE set to median of manufacturing
    theta = theta_nm*[ones(2,1);zeros(K_man,1);ones(K-2-K_man,1)]+[zeros(2,1);B(:,12+j);zeros(K-2-K_man,1)]; %trade elasticities
    gamma = [zeros(2,1);B(:,1+j);zeros(K-2-K_man,1)]; %Heterogeneous scale elasticities, non-manufacturing set to zero
    rho=C(j); %Upper tier elasticity of substitution
  
    for i=1:N
        %Defining data vectors and parameters
        V = V_i(i); %total value added
        X=exports_ik(i,:)'/V; %sectoral exports (divided by value added)
        D= D_i(i)/V; %Aggregate trade balance (divided by value added)
        beta = squeeze(beta_nk(1,i,:)); %sectoral expenditure shares
        lambda= dom_lambda_ik(i,:)';  %own expenditure shares by sector
        V_k = V_ik(i,:)'/V; %value added by sector (divided by value added)
   
    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Optimal trade and industrial policy
        %Specifying taxes and subsidies
        t_k=1./(theta+1);
        s_k=gamma;
        %Solving the model
        model_output=gfop_soe_ces_solver(V_k,D,X,lambda,beta,gamma,theta,rho,t_k,s_k,K);
        gfop_hte_0(i,1)=model_output.welfare;
    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Trade policy only
        %Specifying taxes and subsidies
        t_k=1./(theta+1);
        s_k=zeros(K,1);
        %Solving the model
        model_output=gfop_soe_ces_solver(V_k,D,X,lambda,beta,gamma,theta,rho,t_k,s_k,K);
        VA_TP=model_output.VA_share;
        gftp_hte_0(i,1)=model_output.welfare;

        counter=counter+1
    end

    %%%%%%%%%%%
    %Computing gains from industrial policy
    gfoip_hte_0=gfop_hte_0-gftp_hte_0;
    gfoip_het(j)=mean(gfoip_hte_0);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Common Trade Elasticities
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Using the estimates from the common (median) TE, move the TE up and down
%and the SE inversely

rho=C(11); %Upper tier elasticity of substitution
theta_base=4.32; %median TE
gamma_base = [zeros(2,1);B(:,1+11);zeros(K-2-K_man,1)]; %Heterogeneous scale elasticities, non-manufacturing set to zero
counter=0;
for j=1:13
    %Parameter vectors
    theta=3+.5*(j-1);
    gamma=gamma_base*theta_base/theta;
    for i=1:N
        %Defining data vectors and parameters
        V = V_i(i); %total value added
        X=exports_ik(i,:)'/V; %sectoral exports (divided by value added)
        D= D_i(i)/V; %Aggregate trade balance (divided by value added)
        beta = squeeze(beta_nk(1,i,:)); %sectoral expenditure shares
        lambda= dom_lambda_ik(i,:)';  %own expenditure shares by sector
        V_k = V_ik(i,:)'/V; %value added by sector (divided by value added)
   
    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Optimal trade and industrial policy
        %Specifying taxes and subsidies
        t_k=1./(theta+1);
        s_k=gamma;
        %Solving the model
        model_output=gfop_soe_ces_solver(V_k,D,X,lambda,beta,gamma,theta,rho,t_k,s_k,K);
        gfop_hte_0(i,1)=model_output.welfare;
    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Trade policy only
        %Specifying taxes and subsidies
        t_k=1./(theta+1);
        s_k=zeros(K,1);
        %Solving the model
        model_output=gfop_soe_ces_solver(V_k,D,X,lambda,beta,gamma,theta,rho,t_k,s_k,K);
        VA_TP=model_output.VA_share;
        gftp_hte_0(i,1)=model_output.welfare;

        counter=counter+1
    end

    %%%%%%%%%%%
    %Computing gains from industrial policy
    gfoip_hte_0=gfop_hte_0-gftp_hte_0;
    gfoip_common(j)=mean(gfoip_hte_0);
end

%Putting it all together with the standard devation and theta
theta_common = 3:.5:9;
sd_theta = B(1,24:34);

output_common = [gfoip_common' theta_common']

output_het = [gfoip_het' sd_theta']

csvwrite('gains_alternative_te_common.csv',output_common);

csvwrite('gains_alternative_te_het.csv',output_het);



